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ABSTRACT 


The goal of this study was to investigate the importance of gravity- 
induced free convection in phase-change materials and thereby contribute 
to the understanding of the behavior and performance of phase-change 
thermal control devices. 

Two theoretical models were developed to predict the thermal response 
of the phase-change material to a given hot plate temperature. A two- 
dimensional pure-conduction model was developed to predict the melting 
of the phase-change material when heat transfer was a function of conduc- 
tion. A combined conduction-convection model, also two-dimensional, was 
developed to predict the phase change phenomena when heat transfer was 
a function of conduction and gravity-induced free convection. Both 
models were solved using explicit finite difference approximations on a 
Digital Equipment Corporation, Model PDP-10, digital computer. 

The experimental equipment consisted of a rectangular cell utilizing 
a heating chamber, an expansion chamber, and a test chamber; a sixteen- 
channel multipoint recorder, and a fluid flow system. The recorder 
monitered hot and cold plate temperatures and interior node temperatures 
at two second intervals. 

A comparison of theoretical temperature profiles and experimental 
temperature profiles is presented for six runs at various angles of 
inclination of the test cell with respect to the horizontal direction. A 
detailed discussion of results is presented. As expected, since the pure 
conduction model neglects convection, variation exists between experimental 
results and theoretical results calculated using the pure conduction model 
at angles of inclination other than zero. At an angle of inclination of 
zero degrees good agreement is obtained between experimental data and the 
pure conduction model. Good agreement is obtained between experimental 
data and theoretical temperature profiles using the conduction-convection 
model. The results show that gravity-induced free convection is an impor- 
tant factor in the melting process and can be predicted using the theo- 
retical convection model developed in this study. 
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INTRODUCTION 


Phase change thermal control techniques have received increasing 
attention, (references 1, 23, 24, 25) in the last several years for 
spacecraft thermal design. Because of inherent advantages of simplicity 
and reliability a passive solid-liquid phase change material can be used 
in the walls of spacecraft as packaging around sensitive electronic 
equipment to absorb or release energy to maintain constant temperature of 
the electronic equipment. However, this system is limited by the heat 
rejection or absorption capacity of the material used. 

A previous study^ 2 ) has determined the property requirements of 
phase change materials in order, that they be good thermal control devices. 
The material should be non-toxic, chemically-inert and stable, noncorrosive, 
have small density variations, and have a high latent heat of fusion. The 
material should also melt in the 50- to 150° F range; n-paraffins with an 
even number of carbon atoms are the most widely used materials for this 
purpose. In this study n-octadecane was used. 

j\n earlier study^^ at the Colorado School of Mines dealt with an 
unidimensial melting investigation of a finite paraffin slab. It was 
concluded that the pure conduction model used did not completely solve 
the phase change problem. Therefore, the present study concerns the effects 
of gravity-induced free convection upon the melting phenomena. 

All phase-change experiments, such as ground tests made in high 
gravity fields, must take into account the effect of gravity-induced free 
convection. Either the experiments must be designed to eliminate convec- 
tion or the convection must be mathematically modeled. It is important 
to determine at what gravity level gravi ty-induced free convection may be 
neglected. This will enable designers of phase-change thermal control 
devices for spacecraft ‘to determine whether or not gravi ty-induced free 
convection is an important design factor under low gravity conditions such 
as periods of thrust. 

Other effects, such as electrically-induced convection or magnetically- 
induced convection, may also be important design factors. Since experi- 
ments to study other effects will be made in a high gravity field, the 
effect of gravity must be determined before effects of other forces can be 
.studied completely and modeled accurately. 
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LITERATURE SURVEY 


There has been a large amount of literature published on the subject 
of melting phenomena and gravity-induced free convection. This literature 
survey deals with only a small portion of the published material. One of 
the main references used in this study is the thesis of P. R. Pujado'*). 

In his thesis Mr. Pujado presented a theoretical model for the uni dimen- 
sional melting of a finite paraffin slab. The theoretical model was 
developed using finite difference methods to approximate the solution 
of the partial differential equations governing the physi cal 'system. The 
finite difference approximations were solved on an IBM-Model 360 digital 
computer. The model solved two-phase, uni dimensional heat conduction 
equations with a moving interface and variable thermal properties. Mr. 
Pujado stated that the theoretical model neglected free convection in the 
liquid phase portion of the system and concluded that the errors in his 
results were probably due to the existence of free convection in the cell. 

Earlier. NorthruD Corporation^ 2 ) conducted a similar study and 
obtained results whicn compared very closely with the work Mr. Pujado did. 

Some of the texts which are good theoretical references for heat., 
transfer and fluid flow are Qar s ^ aw anc * 0aeger(3), Rohsenow and Choi' 4 ', 
and Schl i ctinq^ 5 ' . LongwelU 6 ^ was used as the basic tneoretical reference 
Tor developing the boundary layer equations in this study. Dusinberre( ' ) 
was used for development of the interface phase-change equation used in 
the finite difference approximations of the theoretical equations. Bird, 
Stewart, and Lightfoot'®) was used as the reference for free convection 
between infinite parallel plates. Vallentine' 9 ) was used as the basic 
ideal flow reference for the development of the ideal-viscous flow model. 

Grodzka and Fan^ 0 ^ listed various areas of study when attempting to 
solve the problem of free convection in phase change thermal control 
equipment. They stated that free convection might be induced through the 
following forces: gravity, surface tension, electricity, and magnetism. 

The majority of work on free convection effects in liquids and gases 
has been done for infinite plate systems. Models for this type of system 
have been developed by Bodoia and Osterle' 1 '), Dropkin and Globed 2 ), 
Dropkin and Somerscales '* 2 ) , Gebhart^), Kohand Price' *5), anc j Samuels 
and Churchill ' ’ 6 ). 

Wilkes and Churchill^) made a study of temperature profiles in a 
closed rectangular system to determine the effects of gravity induced 
convection. The theoretical model was developed from the basic equations 
of motion, energy, and continuity; a two-dimensional approach precluded 
the study of turbulent flow. The system of equations was solved by an 
implicit alternating-direction technique developed by Peaceman and 
Rachford08). Instabilities in the numerical solutions were noticed above 
certain Grashoff numbers. 
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Other closed cell convection studies have been made by Bellamy-Knights^) 
and Fromm( 2 0). 

Various papers have also been published which deal with the melting 
of finite slabs. Chi-Tien and Yin-Chao YenUU developed approximate 
theoretical solutions for temperature distributions and melting rates when 
the mode of heat transfer was natural convection induced by buoyancy forces. 
Numerical solutions for various ice-water systems were given. Goodman and 
Shea( 22 ' used a series solution to solve the problem of uni dimensional 
melting in a finite slab.. Other works on phase change phenomena include . 
Bannister and Bentilla' 2 ^; Ukanwa, Stermole, and Golden^ 24 '; and Shah' 2 ' 5 '. 

Papers have also been published which discussed Other courses of free 
convection besides gravity-induced convection. Emery ^6) h as studied 
magnetically induced convection. Pearson^ 2 '' and Nield(28) have studied 
the effects of interfacial tension upon convection. Chandrasekan 2 ^) 
studied surface tension effects, rotational effects, and magnetic effects on 
convection patterns. 
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THEORY 


Two separate theoretical models are developed in this study. The 
first model predicts the transient temperature distribution in the system 
when heat transfer is a function of conduction. The. second model predicts 
the transient temperature distribution in the system when heat transfer is 
a combined function of conduction and of gravity-induced free convection. 

The test material used in the research was n-octadecane. The physical 
properties of n-octadecane are given below and in figure 1. The values 
were obtained from Pujado's thesis^ ' , Northrup’s Final Report^), and the 
Data Book on Hydrocarbons . ( 30 ) 

Density 

Solid phase = (-.0008336) T + 1.0918, gm/cc 
Liquid phase = (-.0012505) T + 1.1316, gm/cc 

Heat capacity 

Solid phase = 2.164, watt-sec/gm °K 

Liquid phase = (.008213) T - 0.14237, watt*sec/gm °K 

Conductivity 

Solid phase = (-0. 50054xl0“ 5 ) T + .002914, watts/cm °K 
Liquid phase = (-0.50054x10“^) T + .002914, watts/cm °K 

Melting point = 300.60 °K 

Liquefaction enthalpy = 243.893 watt-sec/gm 

A diagram of the system is given in figure 2. 


Pure Conduction Model 


The equations governing the two-dimensional problem have been developed 
by Carslaw and Jaeger(3). These equations are given below. 


Solid phase 


T = - k S— - 
t Ps c ps 


( T xx + T yy) 


( 1 ) 


Liquid phase 


v 


k a 


h c p 


( T + T ) 
' xx yy 


( 2 ) 




255.2 366.3 477.4 


TEMPERATURE (°K) 
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The boundary conditions for the system under study are 

0 y = 0, T = t d 

@ y = d, T = Tf 

® y = 2y 0> T = T 0 

(3 x = L and x = 0 

Liquid phase 

T = ( T p - T f ) + T f 


Solid phase 

T = y_:.d_ 

2y -d. 

J o 

The initial condition is 

0 6 = 0, T (x,y,0 ) = T. 

The method of excess degrees is used to predict the phase change 
phenomena. When the theoretical solid phase temperature of a volume ele 
ment exceeds the melt temperature a fictitious temperature, (T s -Tf), is 
calculated; when the fictitious temperature, summed over time and multi- 
plied by the heat capacity at the phase change temperature, has a larger 
magnitude than the liquefaction enthalpy, then the volume element has 
changed phase. Figure 3 gives the nodal network diagram. 

An explicit, forward-di fference finite-difference method is used to 
approximate the partial differential equations. The finite difference 
operators are 

T t = T* (n,m) - T(n,m) 

A6 

Txx = T(n+l,m) - 2T(n,m) + T(n-l,m) 

(Ax) 2 

Tyy = T(n,m*-T) - 2T(n,m) + T(n,ro-1) 

( A y )2 


T f - T 0 ) + T 0 




The finite difference operators are substituted into equations (1) 
and (2), and the resulting equations are rearranged to the final form 
given below: 



n+1 , 
m-1 


n+1 ,m 


O 

n+1 , 
m+1 




,m-°2 


T 

Ax 

L 


n,m-l 


n,m+l 


H 


O 

n,m+2 


n-1, 

id-2 


O 

n-1, 

m-1 


O 

n-1 ,m 


n-1, 

m+1 


n-1, 

. _ QL+2 


O 

n-2 , 
m-1 


n-2,m 


r-Ai— i 


o 

n-2, 

m+1 


Figure 3 - Nodal Network Diagram 
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T*(n,m) = T(n,m) (1- 
ks A6 


-2. A6 . k s 


PgCps(Ax) 


T (n+1 ,m) + 


(Ax) PsCp S (Ay) PsCp S 

ks A0 


) 


P sCps(Ax) 


J (n-1 ,m) + 


ks A9 


P S^ps( A ^ 2 

Liquid phase 


, k A6 

T (n jirt+1 ) + -s- 


P s c ps ( A y)‘ 


T(n,m-1) 


(3) 


T*(n,m) = T(n,m) (1- 


2A0k£ 


2Aei<£ 


( A x) 2 p£C p£ (Ay) 2 P£C 


pi 


T(n+l,m) + JlM! 

P £ C f% (A*) 2 p £ C pl ^ Ax ) 

k £ A9 tj. , k $ &e t/. _ n 

J(n,nH-i) + — ^ ; — -J(n,m-i) 

P£Cp£(Ay)^ p £ C p£ G\y) 


(4) 


The following two equations give the stability requirements for the 
finite-difference solution: 


Solid phase 
1 - 


2A0Rs 


2 A0 ks 


> 0 


p s c ds( A x ) 2 p s Cps 


s°ps' 

Liquid phase 

2A0 k£ 2A0 k£ > 0 


1 - 


PjjCp l ( Ax ) 2 P£C p£ (Ay)' 


(5) 


( 6 ) 


The listing of the computer program solving the finite difference 
equations is given in Appendix D. 


Combined Conduction and Convection Model 


Gravi ty-induced free convection affects the heat transfer in the 
liquid phase of the test material. The changes are caused by flow due to 
density variations resulting from the temperature gradient. 

The physical situation under consideration is completely described 
by the following equations with appropriate boundary conditions: momentum 

equation, energy equation, continuity equation, and an equation of state. 

The initial theoretical work was directed toward a numerical solution 
of the above equations and of modified forms of the above equations. 
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However, stability problems were encountered in all numerical solutions. 
Because of these problems an approximate theoretical model has been 
developed. The development of the above equations and a discussion of 
the numerical solutions investigated is given in Appendix B. 

The energy equation for the solid phase is given by equation (1). 
The energy equation for the liquid phase is 

PL = ji L ( t . T ) (?) 

D0 P£C p ^ 1 xx + 1 yy J 

In expanded form the energy equation is 


u T 


v T. 




( T 


PH 


xx 


yy 


(8) 


Substituting in finite difference operators , making the assumption 
that velocity values are those at the node being solved, and rearranging, 
the energy equation becomes 


T*(n,m) = T(n,m) (1 

+ T(n+1 ,m) 

+ T(n-1 ,m) 

+ T(n,m+1) 

+ T(n,m-1 ) 


. Oh A 6 


2kH A6 


' (Ax)* 


u (n ,m) AB 
2 A X 


(JLn.Aex 
Tfi X/ 


+ u (n,m) 
2 Ax 


/ k 9 A8_ 

'TaTF 


V (n,m)AB ) 
2 AY 


( k l A6 
( (Ay) 2 


. v (n,m) A8 ) 
2 Ay 


Ay) 2 


(9) 


The stability criterion from the solid-phase energy equation is 

1 _ 2A9k s _ 2 A9 k s > 0 (10) 

P 5 C ps (ax ) 2 psCps ( ay ) 2 

The stability criterion from the liquid-phase energy equation 



2 AG k 


P C 

TP* 



1 


> o 


(ID 
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u:. (n,m) 
v . (n ,m) 


_2ko > 

"tT ? 1 

2k p > 

P£ C n Ay 
Pi 


0 

0 


( 12 ) 

(13) 


Equations (12) and (13) give maximum stable values of velocity that 
can be used in the finite difference solution of the liquid-phase energy 
equation. Taole I gives tabulated stability criteria for various values 
of Ax and Ay. 


Table I 

Stability Criteria 


(ax) 

CM 

(Ay) 

CM 

Mi 

Sec 

A 0 s 

Sec 

u 

CM/Sec 

V 

CM/Sec 

.3175 

.15875 

12.57 

12.25 

.00570 

.01112 

.2539 

.2539 

18.26 

18.17 

.00695 

.00695 

.2539 

.15875 

10.263 

10.00 

.00695 

.0112 

.15875 

.079375 

2.855 

2.781 

.01112 

.02224 

.079375 

.079375 

1.78 

' 1.74 

.02224 

.02224 


The velocity profile used in the convection model is an approximate 
profile obtained by combining an ideal flow solution for flow in a 
cul-de-sac region^) with a viscous flow solution for flow between 
infinite parallel plates. A maximum velocity is imposed on the ideal- 
viscous flow pattern, using either a driving velocity calculated from the 
buoyancy force term or the maximum velocity calculated from the liquid- 
phase energy equation stability criteria. 

The velocity profile for unidimensional flow between infinite vertical 
plates is developed by Bird, Stewart, and Lightfoot. (8) The velocity pro- 
file i s 

y =PJL gkjbl (7? 3 -??) (14) 


where rj = 



12 


Since a maximum velocity has been determined the velocity profile 
should be given as a function of the maximum velocity. The equation 
for velocity now becomes 


v= 


v max 



(15) 


where 



(from centerline) 


The ideal flow velocity profiles are developed by the use of complex 
variable transformations; the flow pattern under consideration is shown 
in the following diagram 


B A 



The basic assumptions made in the conformal transformation process 
are that (1) the flow pattern being studied is irrotational flow of a 
perfect fluid and that (2) the complicated flow pattern can be trans- 
formed by use of complex variables into parallel, uniform flow. 

The flow pattern shown above is assumed to be a complex z-plane 
flow within straight-wall boundaries. Since we are investigating ideal 
flow with a simple polygon, a Schwarz-Christoffel transformation™) may 
be used to obtain a parallel uniform flow pattern. 

If the polygon is in the z-plane and the new plane is the t-plane, 

then 


|| > A’ta-tr^Cb-t)^ (16) 

where P? = complex constant 

a,b,c = real constants in ascending order of magnitude 
external deflection angles of the polygon 
for the flow pattern under study 

-VM- -e/7r= -Un- - | 
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The boundary conditions for the transformation are 

@A , t = - oo 
@B, t = - ! 

@C, t = + 1 
@D, t = +«> 

Therefore 

r dt 

z = A 1 J V(l-t) (-1-t) +B 1 


(17) 


or z = a! cosh t + 


(18) 


Applying the boundary conditions 
@c t = 1 , z = o 
o = A^ cosh (1 ) +b! 

0 = B 1 

@B , t = -1 , z = il 
i i = A 1 cosh“^ (-1) 

1 i = A^i7r 
Al = iln 

Therefore 

z = i/7i cosh (t) (19) 

According to definition as given by Vallentine (^) parallel, uniform 
ideal flow should have a source at - 00 and a sink at . However, the 
above t-plane flow pattern gives a source at + 00 and sink at - °° . 
Therefore w = -t, where w is a new complex plane, and 

w = - cosh in z. ) 

l 

, where z = x+iy 

U =-cosh (ZL* iJLZ) 

i i 

or w = r cosh(7rx)cos(7ry) - i sinh(7tx) sin (tty) 

1 1 1 1 


(20) 
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By definition of complex flow 


9 = -coshroc cos7ry (21) 

i i 

ip = -sinh 7r x sin n y (22) 

i i 

The stream function is defined by the following equations 
u = ip y (23) 

v = -^ x (24) 


Substituting into equation (22) and differentiating we obtain 

u = -2£sinh ( — x) cos ( — y) (25) 

i i i 

and v = —cosh ( — x) sin ( — y) (26) 

i i i 

The liquid phase is split into three regions; as shown below. 


Region 

Region 

Region 

I 

II 

HI 

IDEAL FLOW 

VISCOUS FLOW 

IDEAL FLOW 


Region I flow is governed by equations (22), (25), and (26) with 
appropriate boundary conditions. Region II is governed by equation 
(15) with a given maximum velocity. Region III is governed by another set 
of ideal flow equations; but assuming symmetrical flow it is not necessary 
to develop the equations for this region. 

The ideal flow regions are coupled to the viscous flow region by 
assuming that the velocities in the viscous flow region are the boundary 
velocities for the ideal flow reqion. All y velocities are zero at this 
boundary. By use of equation (22) values of the stream function may be 
calculated at the ideal flow-viscous flow boundary. Since there can be 
no flow across a line of constant stream function, velocities in the ideal 
flow region can be related to boundary velocities at calculated values 
of the stream function at the boundary. By this method a pseudo-viscous 
flow pattern can be imposed on the ideal flow regions. 

The actual liquid phase in the test cell does not have a constant 
depth. The velocities calculated by the ideal flow-viscous flow model 
are imposed on the liquid phase at a point by assuming the depth of the 
liquid phase at that point to be the depth of the cell in the ideal 
flow-viscous flow model. 

The flow pattern calculated is only an approximation, but with the 
small magnitude of allowable velocities calculated by stability criteria 
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the velocities calculated should give fairly accurate flow patterns in 
the liquid phase. 

A listing of the conduction-convection model computer program is 
given in Appendix E. 

An additional source of error exists in the finite difference solution 
of equation (8). The finite difference approximation neglects the second 
order partial differential with respect to time, and considers it be a trunca- 
tion error to be included with higher order truncation errors. 

If the second order partial is kept equation (8) now becomes 


V 


At 


T tt + U T x + v T y -ft < T xx +T yy > 


(27) 


When equation (8) is rearranged to solve for T and then substituted 
into equation (27) the resulting equation is 2 

t 4-viT + ttt = t k At u \ j . / k _ \ y 

t T x T y ~ pc p " 2 T xx pc " 2 T yy 


- At uv T, 


yx 


(28) 


The magnitude of the numerical dispersion coefficient is approximately 
one-half the value of the thermal diffusivity coefficient for values of physi- 
cal properties .time step, and velocity levels used in the numerical solution. 
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EQUIPMENT AND PROCEDURE 


The test cell, figure 4, consisted of a rectangular test chamber, 
a heating chamber, and an expansion chamber. The test chamber was 
2.54 cm by 12.7 cm by 12.7 cm. Fourteen i ron-constantan thermocouples 
were placed in the test chamber for the purpose of recording temperatures. 

A valve was placed in the heating chamber, which connected the test 
chamber and expansion chamber. The expansion chamber was included in 
the cell design because the cell was originally designed for use in a 
vacuum; this fact necessitated the use of a completely closed system. 

The heating chamber consisted of a rectangular copper cell with an entrance 
port on each of the vertical walls and two exit ports through the top 
plate of the chamber. 

Although the cell was designed and built as a completely closed 
system, with all edges sealed, leaks were encountered when the cell was 
first used. The problem was overcome by the use of an epoxy coating on 
all joints, and by making experimental runs at lower temperatures than 
originally planned. 

The heating system, figure 5, consisted of a constant temperature 
bath, a centrifugal pump, 6 tygon lines connecting the bath, pump, and 
test cell. The constant temperature bath was a 5-liter glass tank, two 
immersion heaters, and a stirrer. The heating fluid used was water. 

The recorder used was a Bristol multipoint recorder. Sixteen channels 
were used, with a two-second interval between points recorded. The ac- 
curacy of the recorder was ±0.417°K. The leads from the test cell were 
connected directly into the recorder. 

The cell filling apparatus, figure 6, consisted of a cell filler 
and a constant temperature bath. The test material was heated by running 
water in coils from the constant temperature bath through the cell filler. 
The bath was kept at 305°K; this temperature was used because a small 
solid-liquid density change in the test material was desired when filling 
the cell. The test material was degassed by the use of a magnetic stirrer. 

Experimental runs were made using the following procedure: 

1. The constant temperature bath was allowed to heat to approximately 
2.22°K higher than the desired hot plate temperature; this procedure 
allowed for the small cooling effect of the cold water present in the 
expansion chamber of the test cell. 

2. The recorder was allowed to run during the heating period in 
order to check the initial steady state temperatures and to monitor the 
heating tank temperature. 

3. When the tank was at the desired temperature the run was initiated 
by turning on the pump. 



Note: All dimensions are in cm 


1 ? 



Figure 4 - Pseudosection of Test Cell: Scale 1cm: 1cm 




Figure 5 - Equipment Flowchart 
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4. The duration of the run was approximately 40 minutes. 

Runs were made at angles of inclination of 0-, 30- , and 60- degrees. 
Runs were not made at higher angles of inclination because the phase- 
change material melted completely in the top part of the test cell. 
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The results from six data runs have been compared with the 
theoretical temperature profiles predicted by the pure conduc- 
tion model program and the conduction-convection model program. 

The results obtained show the effect that gravity has upon the 
phase change process . The results show that the pure conduction 
model cannot predict actual temperature profiles in the phase 
change material if gravity-induced free convection is present. 
However, good agreement is obtained for results comparing ex- 
perimental data and the conduction- convection model. 

Figure 7 prgsents results for y = 0.635 cm at an angle of 
inclination of 0°. The spread in experimental results is due to 
the presence of air bubbles in the liquid phase. At an angle of 
inclination of 0 the pure conduction model predicts the experi- 
mental temperature profiles closely, considering that the air 
bubbles do affect the experimental temperature profiles . 

Figure 8 shgws results for y = 1.27 cm at an angle of 
inclination of 0°. Good agreement is obtained between experi- 
mental data and the pure conduction model temperature profile. 

Figure 9 presents a comparison of experimental data to the pure 
conduction temperature profile for the same run at y = 1.905. 

Again good agreement is obtained between the theoretical tem- 
perature profile and experimentally measured temperatures. 

Figures 10, 11, and 12 present comparions of experimental 
data to theoretical pure conduction temperature profiles for 
Run 7 at y equal 0.635 cm, 1.27 cm, and 1.905 cm, respectively. 

The angle of inclination is again 0°. The 0.635 cm results 
again show spread. The main cause of thy spread is due to air 
bubbles . In figure 7 and figure 10 the thermocouple at the 
middle of the cell records a lower temperature than the thermo- 
couples at the ends of the cell. Part of this difference may 
be due to the fact that the plexiglas has a slightly higher 
thermal diffusivity than the n-octadecane . Therefore, the 
results are also showing end effects due to higher temperatures 
in the plexiglas walls. The 1.27 cm and 1.905 cm results again 
present good agreement between experimental results and 
theoretical prediction. 

Figures 13, 14, and 15 present a comparison of experimental 
data to pure conduction temperature profiles for Run 3 at y 
equal 0.3175 cm, 0.635 cm, and 1.27 cm, respectively. These re- 
sults are presented to show qualitatively the effects of free 
convection. In figure 13, the thermocouple at x = 2.8575 cm 
takes the longest time to melt; the thermocouple at x = 11.1125 cm 
takes the shortest time to melt. The thermocouple at x = 6.6675 cm 
melts at an intermediate time. This trend is to be expected 
when gravity- induced free convection is present. The depth of 
liquid should be larger at the top of the cell than at the bottom; 
as measured along the long axis of the cell the large x-positions 
will be referred to as the top of the cell, intermediate x-positions 
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as the middle of the cell, and small x-positions as the bottom 
of the cell. Since the pure conduction model predicts the 
same temperature profile for all three x-positions it is 
evident that the model cannot predict the experimental temperature 
profiles obtained. 

The same trend is present in figure 14; but in this case 
only the thermocouple at x = 10.16 cm deviates from the others 
during the duration of the experimental run. The pure conduc- 
tion model cannot predict this deviation in experimental tem- 
peratures . 

Figure 15 again shows the deviation of experimental results 
from the pure conduction model prediction. At this y position 
the time required for the deviation in temperature to occur is 
longer than in figures 13 and 14; but still it is evident that 
the pure conduction model does not predict the temperature pro- 
files when convection is present. 

Figures 16 and 17 present theoretical and experimental com- 
parisons for Run 4 at y = 0.3175 cm and y = 0.635 cm. These 
results show the same trends as figures 13 and 14. They are 
presented to show that convection occurs in more than just one 
run. The pure conduction model does not predict the effect of 
convection evident in experimental data. 

The results presented in figures 7 through 17 show that 
the pure conduction model will predict experimental temperature 
profiles when convection is not present, but will not predict 
experimental temperature profiles when convection is present. 

Figures 18 through 44 present a comparison of experimental 
data to both the convection model and pure conduction model for 
angles of inclination of 30 and 60 . If only one theoretical 
curve is shown in a figure it indicates that both theoretical 
models predict approximately the same temperature profile. 

Figures 18 through 26 show results for Runs 8 and 10, 
both made at an angle of inclination 60 . Some spread in ex- 
perimental results is evident; the spread is due to the presence 
of air bubbles. Figures 18, 19, and 20 give results at small 
x-positions. Figure 19 shows that the convection model predicts 
the experimental results more closely than the pure conduction 
model. It should also be noticed that the convection model 
predicts a lower temperaturyprof ile than the pure conduction 
model. Figures 20 and 21 show good agreement between both 
theoretical models and experimental data. 

Figures 21, 22, and 23 present results at intermediate 
x-positions . Both models predict approximately the same tem- 
perature profile; there is a small difference in figure 21, 
due to different nodal sizes used in the models. Both models 
predict the experimental results, taking into account some 
spread in experimental results, due to air bubbles (see 
Figure 21) . 
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Figures 24, 25, and 26, at the larger x-positions, show 
the largest deviations between the pure conduction model and 
the convection model. In all three figures the convection 
model has better agreement with experimental data than does 
the pure conduction model. Deviation between figure 19 shows 
that the convection model predicts the experimental results 
more closely than the pure conduction model. It should also 
be noticed that the convection model predicts a lower tempera- 
ture profile than the pure conduction model. Figures 20 and 
21 show good agreement between both theoretical models and 
experimental data. 

Figures 21, 22, and 23 present results at intermediate 
x-positions . Both models predict approximately the same 
temperature profile; there is a small difference in figure 21, 
due to different nodal sizes used in the models. Both models 
predict the experimental results , taking into account some 
spread in experimental results, due to air bubbles (see figure 
21 ) . 


Figures 24, 25, and 26, at the larger x-positions, show 
the largest deviations between the pure conduction model and 
the convection model. In all three figures the convection 
model has better agreement with experimental data than does 
the pure conduction model. Deviation between experimental 
data and the convection model is caused by the reduced gravity 
field used in the convection model and imposed by velocity 
stability criteria in the energy equation for the model At 
this end of the cell the convection model predicts higher 
temperatures than predicted by the pure conduction model. This 
trend agrees with experimental results obtained. 

Figures 27 through 35 present results obtained for a 30° 
run. Figures 27, 28, and 29 represent thermocouples located 
at small x-positions. The y = 1.27 cm and y = 1.905 cm 
theoretical curve, both convection and conduction models, 
agrees very closely with the experimental results. At 
y = 0.635 cm the convection model agrees closely with experi- 
mental data; the pure conduction model deviates from the 
experimental data and predicts higher temperatures than obtained 
experimentally. Figures 30, 31, and 32 present a comparison 
of theoretical and experimental results at intermediate x- 
positions . At all y positions good agreement is obtained 
between theory and data; both models predict the same tempera- 
ture profiles at all y positions. Figures 33, 34, and 35 
represent thermocouples located at large x-positions. The 
convection predicts all temperature profiles closely; while 
deviation occurs between the pure conduction model and experi- 
mental data at y = 0.635 cm and y = 1.27 cm. The convection 
model predicts higher temperatures than does the pure conduction 
model. 



24 


Figurgs 36 through 43 present results obtained for 
another 30 angle of inclination run. The results from this 
run are presented separately because the initial temperature 
is different from run 4. Figure 36, 37, and 38 show a compari- 
son of experimental data and theoretical temperature profiles 
for small x-positons . As in the previous cases the y = 1.27 cm 
and y = 1.905 cm results show good agreement between theory 
and data. At y = 0.635 cm the convection model predicts a 
lower temperature profile than does the pure conduction model; 
again the convection model temperature profile is closer 
to the experimental data than the conduction model temperature 
profile. Figures 39, 40, and 41 present a comparison of 
theoretical and experimental results at intermediate x-positions. 
The theoretical results, both conduction and convection, agree 
very closely with experimental data. There is some deviation, 
at y = 0.635 cm, due to air bubbles. 

Figures 42, 43, and 44 show comparisons of theory and 
data for large x-positions. At all y-positions the convection 
model more closely predicts the experimental data than does 
the pure conduction model. At y = 0.635 cm and at y = 1.905 cm 
the experimental phase change takes place sooner than theoret- 
ically predicted by the convection model; the deviation is due 
to stability limitations and numerical dispersion effects. 

Figures 45, 46, and 47 present theoretical results show- 
ing the effect of velocity level upon the shape of the solid- 
liquid interface. At v m x = 0.0107 cm/sec the maximum difference 
in interface position between the top, x = 10.16 cm, and bottom, 
x = 2.54 cm, is 0.762 cm. At v = 0.00535 cm the maximum 
difference is 0.278 cm. At v ma = 0.002675 cm the maximum 
difference is 0.056 cm; this Sitference is negligible since the 
y used in the calculation was 0 .076 cm. At this theoretical 
velocity level the pure conduction model can be used to predict 
temperature profiles and gravity level. 

Figure 48 is included to show the effect of velocity level 
upon temperature profiles. At x = 10.16 cm, y = 0.635 cm the 
pure conduction model has the largest deviation from a typical 
experimental temperature profile. As the velocity level in- 
creases the convection effect becomes larger and the theoretical 
profiles more closely predict the experimental profile. Since, 
due to stability criteria, 0.0107 cm/sec is the largest allow- 
able maximum velocity the experimental profile cannot be exactly 
predicted. There is no effect due to convection at the inter- 
mediate x-position, figure 48-b. At the low x-position, figure 
48-c, the results show that the convection model will predict the 
experimental results closely even at low velocity levels. 

Figure 49 presents results from theoretical computer 
runs made to determine numerical dispersion effects. 

At the middle of the test cell, figure 49-b, there is no 
apparent effect of numerical dispersion. At the top and 
bottom of the cell, figures 49-a and 49-c, the effect shows 
the same trends as varying the velocity level, but with the 
time step range used in this study the effect is smaller than 
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varying the velocity. The solutions are not convergent 
with the time steps used, but due to computer limitations 
smaller time steps could not be used. 

No comparisons between experimental data and the convection 
model is presented for y = 0.3175 cm in this discussion of 
results. Because of air bubbles large deviations occur between 
data and theory. These results are presented in Appendix C. 



FIGURE 7 - Comparison of data to pure conduction temperature 

profile for Run 6 at y = 0.635 cm 
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FIGURE 8 - Comparison of data to pure conduction temperature 

304.0 - 

profile for Run 6 at y = 1.27 cm 
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FIGURE 10 - Comparison of data to pure conduction 
temperature profile for Run 7 at y = 0.635 cm 
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Figure 17 - Comparison of data to pure conduction temperature profile 
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Figure 21 - Data versus theory for y = O .635 cm, x = 5.715 cm 

0 . = 60°, O = Run 8 , 0 = Run 10, = theory, convection 

= theory, pure conduction 
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Figure 22 - Data versus theory for y = 1.27 cm, x = 6.0375 cm 

a = 60°, O = Run 8 , Q = Run 10, = theory, convection 

= theory, pure conduction 
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Figure 23 - Data versus theory for y = 1.905 cm, x = 6.35 cm 

a = 60°, Q = Run 8 , 0 = Run 10, = theory, convection 

= theory, pure conduction 
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Figure 24 - Data versus theory for y = O .635 cm, x = 10^_ljS- 
d = 60°, O = Run 8, ^ ^ 

□ = Run 10 © Q 


-= theory, convection 
•= theory, pure conduction 


12 

00 

TIME 

( SEC ) 


Figure 25 - Data versus theory for y 
a = 60°, © = Run 8, Q = Run 10 

= theory, convection 

= theory, pure conduction 
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Figure 33 - Data versus theory 

for y = 0.635 cm, x = 10.l6jyp. 

a = 30°, O = Run 4 
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Figure 34 - Data versus theory for y = 1.27 cm, x 
a = 30°, O = Run 4 
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Figure 35 - Data versus theory for y = 1.095 cm, x = 10.795 cm 
a = 30°, O = Run 4 
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Figure 40 - Data versus theory for y = 1.27cm, x = 6.0325cm 
a = 30°, O = Hun 3 


•= theory, both models 



Figure 41 - Data versus theory for y = 1.905cm, x = 6.35cm 
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Figure 48 - Effect of velocity level on 


profile 
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at 

y = 0.635 
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Legend : 
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Figure ^9 - Effect of numerical dispersion on 
temperature profile 

(a) at y = 0.635 cm, x = 10.160 cm, a = 60 

(b) at y = 0.635 cm, x = 5*715 cm, a = 60 ' 

(c) at y = 0.635 cm, x = 1.905 cm, a- 60 


Legend: 
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CONCLUSIONS 

The following conclusions were drawn from this study: 

1. The study has demonstrated that gravity-induced free con- 
vection greatly alters the melting interface profile and 
the temperature profiles of individual nodes within the 
phase change material. 

2 . The results show that the pure conduction model cannot pre- 
dict gravity effects. In cases where convection is not 
present the pure conduction model does a good job of model- 
ing the phase change process. 

3. The convection model shows good agreement with experimental 
data. The velocity profile model is accurate enough to 
give temperature profiles to be used in the preliminary 
design of phase- change thermal control packages. 

4. It was not possible to determine an approximate gravity 

level. Methods have been suggested^^ that relate a 
critical Rayleigh number to the velocity level. However, 
literature values for the critical Rayleigh number vary 

(31) 

from 600 to 4000 for conditions similar to the one 

being studied. Therefore, any Rayleigh number in the range 
could be used to determine a gravity level which would 
justify the theoretical model. However, there is no basis 
for choosing a given Rayleigh number. Because of the un- 
certainty in the value of the Rayleigh number, we cannot 
state that the results correlate with a known gravity 
level. We can, however, evaluate convection effects in 
terms of the maximum velocity parameter and thus investigate 
convection effects in phase change devices. 

5. Within the accuracy of numerical convection solutions ob- 
tained, the results indicate that gravity-induced free 
convection may be neglected at velocity levels less than 
0.002675 cm/sec. 

6. The experimental investigation has shown that air bubbles 
have a large effect upon the performance characteristics 
of the phase change material. 

7. Sources of error in the numerical solution are: 

(a) numerical dispersion effects 

(b) arbitary flow profiles 

(c) stability limitations 

(d) constant maximum velocity 

8. The method of solution is an initial solution to the 
phenomena of gravity-induced free convection in the 
liquefaction of a phase change material. The study has 
indicated problem areas in numerical solutions of the 
physical situation and has shown that further work is 
needed if complete solutions to the problem are to be 
realized. 
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RECOMMENDATIONS 

The following recommendations are made as a result of this 

s tudy . 

1. Because of the large effect that air bubbles have upon 
the performance of the phase change material, further 
investigation should be made on methods of proper degass- 
ing of phase change materials. 

2. A theoretical investigation should be made to determine 
the proper form of finite difference approximation needed 
to introduce an equation of state for density into the 
mathematical model for velocity. 

3. An experimental study using tracer materials in the phase 
change system should be undertaken to determine the actual 
shape of convection induced velocity profiles in liquifaction 
and solidification phenomena. The measurement equipment 
should be photographic or microphotographic equipment. 

4. A study should be made to determine the magnitude of gravity- 
induced free convection effects upon PCM performance in 

a PCM-filler system. 

5. Because of the problems encountered in the theoretical 
modeling of the convection phenomena, a material investigation 
study should be made to determine whether or not low density 
polymers, which demonstrate a solid-solid phase change of 
proper magnitude in the proper temperature range , would make 
feasible phase-change materials. This type of material could 
be modeled using only a pure conduction model. 

6. Further studies should be made, using a phase change material 
which has physical properties such that the stability criteria 
developed in the finite difference approximation for the 
vorticity equation is the governing stability criteria 

for the theoretical model. Then, given a computer with 
large enough memory capacity, the convection-conduction 
system could be accurately modeled. 
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Nomenclature 

Definition: 

Given that s = f(x,y,t), then the following definit- 
ions are true: 

_6_s _6_s _6_s _ 6^s _ 6^s 

s t 6t» S x 6x» S y 6 y ’ s tt 6t 2» S xx " 6x 2> 

e£s 6 2 s 

yy “ 6y 2’ *y - 6x6y ■ 

Parallel Plow Model: 

b = one-half of distance between parallel walls, cm 

p 

g = acceleration of gravity, cm /sec 
AT = temperature gradient between parallel walls, °K 
y = distance from centerline, cm 

V = y/b 

0 = coefficient of thermal expansion, °K~* 
p = density, grams/cubic cm 
p = viscosity, gm cm”*sec“* 

Ideal Flow Model: 

a,b,c = real constants in ascending order of magnitude 
A',B' = complex constants 
t,w,z = complex planes 

u = x-direction velocity, cm/sec 
v = y-direction velocity, cm/sec 
x = spatial dimension in complex z plane, cm 

y = spatial dimension in complex z plane, cm 

— 1 

9 = potential function, sec 
fjf = stream function, sec -1 
Finite Difference Models: 

c = heat capacity, watt® sec gm"* °K“* 

■t' __ -1 

AH = latent heat of fusion, watts sec gm" 
k = thermal conductivity, watts cm"* °K~* 

T = temperature, °K 
u = x-direction velocity, cm/sec 
v = y-direction velocity, cm/sec 
x = spatial direction, cm 
y = spatial direction, cm 
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a = angle of inclination, degrees 
p = density, gm/cc 
B = time, sec 
A6 = time step, sec 
v = kinematic viscosity, centistokes 
Subscripts: 

f,i,l,o,p,s = fusion, initial, liquid, cold plate, 
hot plate, solid, respectively 


Superscript: 

* = at new time step 
Pure Conduction Computer Program: 

AHP = latent heat of fusion, Btu/lb 

AKL = liquid phase thermal conductivity, Btu (ft sec 0 F)~' 1 ' 
AKS = solid phase thermal conductivity, Btu (ft sec 0 F) -1 


(lb °F) _1 


CPL = liquid phase heat capacity, Btu (lb 0 F) -1 
CPS = solid phase heat capacity, Btu 
DT = time step, sec 

DX = spatial increment, x-direction, inches 
DY = spatial increment, y-direction, inches 
RHOL = liquid phase density, lb/ cubic ft 
RHOS =. solid phase density, lb/ cubic ft 
T = temperature, °F 
TAU = initial time, sec 


TE 

TF 

TO 

TP 


= excess degrees, F 

= fusion temperature, °F 

= cold plate temperature, °F 

o T 


hot plate temperature, ~F 
Combined Model Computer Program: 

AHF = latent heat of fusion, Btu/lb 

2 

AL = liquid phase thermal diffusivity, ft /sec 
AS = solid phase thermal diffusivity, ft /sec 
DX = spatial increment, x-direction, ft 
DY = spatial increment, y-direction, ft 
R = phase indicator 
T = temperature at. old time step, °F 
TD = time increment, sec 
TI = elapsed experimental time, sec 
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TF = phase change temperature, °F 
TN = temperature at new time step, °F 
TO = cold plate temperature, °F 
TP = hot plate temperature, °F 
UX = x-direction velocity, ft/sec 
UY = y-direction velocity, ft/sec 
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APPENDIX A 

The experimental data for this study is given in the 
following study: 

Bain, R. L., "The Effect of Gravity-Induced Free 
Convection Upon the Melting Phenomena of a Finite 
Paraffin Slab for Thermal Control," Thesis No. 

T 1319, Colorado School of Mines, Golden, Colorado, 
1971. 
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Discussion of Other Theoretical Models 


Other theoretical models were developed in this 
study, but due to stability requirements and other numer- 
ical difficulties it is not possible to use these models 
as solutions at the present time. Only one development 
is presented in this section; this model seems to have 
the best possible application to the problem of gravity- 
induced free convection effects in solid-liquid phase 
change. 


The basic equations used in this development are the 

/ O \ 

equations of motion, energy, and continuity. ' These 
equations are given below: 

Motion 

u t +uu x +vu y = v ( u xx +u yy )-q(T— T a ) (29) 

where 

q=- agsina 

^a 

and, 

v t +uv x +vv y = ^(v xx +v yy )-r(T-T^) ( 30 ) 

where 



Continuity 


V T y = 0 


Energy 


T . +uT +vT 
t x y 


(3D 

(32) 


The stream function and vorticity, defined below, 
are introduced after the equations of motion are differ- 
entiated, the u equation with respect to y and the v 
equation with respect to x, and the equation of the 
difference between the two developed. 

Stream function 

u ■ % 

V = -P x 
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Vorticity 

X. = - 9 2 = -u + v (3M 

J A 

The resulting vorticity equation is given below: 

h + <^>x + < rt >, “ a xx + \y> - « T y ' rT x <35) 

The equation of continuity is inherently satisfied 

by the vorticity equation. Equations 32, 33, and 35 are 

then put into finite difference form; these equations 

are used, with appropriate boundary conditions, to model 

the liquid phase. Methods of finite difference formulation 

( 17 ) ( 32 ) 

of these equations are presented by Wilkes and Fromm. ' 

There are two problems present in the theoretical 
approach which make the solution meaningless at the present 
time. First, the base temperatures, T q and T^, in the 
gravity approximation are not well defined. Average 
values of the liquid phase temperature give temperature 
driving forces which cause velocities to exceed stability 
criteria after one time step. Therefore, before any 
finite difference solution using this gravity approximat- 
ion can be made further investigation needs to be made 
in two areas. The first area to be studied is the correct 
determination of values of T and T' to be used in the 

3 3 

gravity approximation. The second area to be studied is 
alternate finite difference formulations of the gravity 
approximation in order to reduce the magnitude of the 
gravity effect in the finite difference formulation. 

An implicit finite difference technique can be used 
to eliminate stability requirements with respect to the 
magnitude of the time step; but a stability requirement 
still exists with respect to the maximum allowable velocity. 
The stability requirements for velocity in the energy 
equations are given by equations 12 and 13". The stability 
requirements for velocity in the vorticity equation are 

u(n,m) - 2v/Ax > 0 


(36) 
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v(n,m) - 2 v/Ayi > 0 (37) 

In the finite difference solution velocities velo- 
cities are calculated from the vorticity and stream funct- 
ion equations, then used in the solution of the energy 
equation. Therefore, to ensure a stable solution the 
properties of the test material must be such that 

v > k/pc p (38) 

Otherwise, velocities may be calculated which are stable 
in the vorticity equation, but which cause the energy 
equation to be unstable. N-octadecane does not exhibit 
this property. 
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APPENDIX G 
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TIME ( SEC ) 



TEMP (°K) 
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TEMP 
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TEMP. (°K) _ TEMP ( K) 
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TIME ( SEC ) 
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APPENDIX D 



PURE CONDUCTION COMPUTER PROGRAM 
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This program was written in FORTRAN-II to solve a two- 
dimensional pure- conduction liquefaction problem. The program 
was run on a CDC-8090 computer. Consult nomenclature section 
for definition of terms. 


Input file: 


1. First card - DX, DY, TD, TAU, AHF , TF , TI 
where the format is (7F10.5), and 


DX [=] 
DY [=] 
TD [=] 
TAU [=] 
AHF [=] 
TF [=] 
TI [=] 


inches 

inches 

minutes 

mi nutes 

BTU/LB 

°F , 

Min 


2. Second card - AT, BT , CT , DT , ET, FT 

where the format is (5E16.8) , and all variables are 
unitless. This is a dummy input file, to be used if 
a polynomial fit of hot plate temperature is desired; 
statement no 12 in the listing will have to be changed 
to polynomial form to use a polynomial fit. The poly- 
nomial is of the form 

T = AT * TAU 5 + BT* TAU 4 + CT * TAU 3 + DT * TAU 2 
P 

+ ET * TAU + FT 

3. Third card - ACT, BCT , CCT , DCT , ECT , FCT 

where the format is (5E16.8), and all variables are 
unitless. The polynomial is of the form 

T = ACT * TAU 5 + BCT * TAU 4 + CCT * TAU 3 + DCT * TAU 2 
o 

+ ECT * TAU + FCT 

4. Sample Input 

Card 1 - 0.250, 0.0625, 0.0, 104.90, 81.50, 0.98 
Card 2 - 0.1E-20, 0.1E-20, 0.1E-20, 0.1E-20, 0.1E-20, 

0.1E-20 

Card 3 - -.109E-05, .11859E-03, -.4566E-02, .6941E-01, 
- . 1157E + 00, . 74487E + 02 


Output file: 

This program will print 

1. Input variables 

2. Time, t, at one minute intervals 

3. Nodal Indicators for all nodes at time t 

Solid = +1 

Solid (at phase change temp.) = +3 
Liquid = -1 
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4. 


Temperatures 


for all nodes at time t. 


Note : 


Subroutines TS1, TS2, and TLl calculate the physical 
properties of n-octadecane as functions of temperature for 
solid, interface, and liquid phase nodes, respectively. For 
a different test material, these subroutines will have to be 
modified. The units used for physical properties in the 
subroutines are given below 

C £=] Btu lb _1 °F -1 
P 

k [=] Btu (hr ft °F) ^ 
p [=] lb (cu. ft.) ^ 



dimension T(22, is > ,cps( 22 , 18 > ,cpk 22 , is > ,aks< 22 , 18 ).,akl< 22,18 ) ,khq 
lS(22,18),RH0L(22,16),R(22fl8},TE(22,l8) 

COMMON TO,DX,DY 7 o 

C READ DATA 

Read i,dx,dy»td,tau»ahf,tf,ti 

1 FORMAT ( 7F 10 , 5 ) 

PRINT 1,OX,OY,TDiTAU,AHF,TF,TI 

DX=DX/12,0 

OYsDY/12,0 

READ 2,AT*BT,CT,DT»ET,FT 
READ 2, ACT^CT.CCT.OCTfECTjFCT 

2 FORMAT (5E16.8) 

PRINT 2, AT,BT,CT,DT.ET,FT 
PRINT 2, ACT,BCT,CCT,DCT,ECT,FCT 
' C SET NODAL INDICATORS AND INITIAL CONDITIONS 
RMAXS0.0 
DO 8 M«1,18 
DO 8 N*l,22 
RIN.MJsl.a 
T ( N , H ) s F C T 
TE(N,M)=0,0 . 

3 CONTINUE 
C CALCULATIONS 

10 TAU=TAU*TO 

12 TP = 120 , 0 

13 TO= ( < < ( ACT»TAU*BCT)*TAU*CCT)*TAU*DCT)«TAU*ECT)*TAU*FCT 
DO 25 N*1 , 21 

T ( N , 1 ) a T P 
25 T { N, 17 ) =TO 
24 FORMAT (2F10.5) 

IF ( 7P*TF ) 29,29,27 

27 DO 26 M?2,16 
WsM 

TU,M)sTF«(W/l6, )*(TF»TO> 

28 T(21#M)*T(i,M) 

DO TO 99 

29 DO 39 M?2 , 16 
WsM 

T (1,M)=TP»(W/16, )»(TP«TO> 

39 TC21,M)=T(1,M) 

99 DO 50 M»2 , 16 
DO 40 N ? 2 , 2 0 
IF ( R ( N , M ) ) 30,100,15 

15 IF <R(N,M>*2,0> 31,100,35 
IF (T(N,M)*TF) 40,16,16 

31 CALL TS1(T,N,M,aSi,AS2,AS3,BSi,BS2,BS3> 

GO TO 36 

35 CALL TS2CT,N,M,AS1,AS2,AS3,8Si,BS2,BS3> 

36 T(N,M)=T(N,M)#(1,0»AS1«2,0*BSi#2,0)*AS2*T(N*1,M)*AS3 # T(Nb1,M)*BS2 # 
IT (N,M*1)+BS3*T ( N , M * 1 ) 

IF(T(N,M>»TF) 40,16,16 

16 SUM=T(N,M)^TF 
SUM=SUM+TE(N,M> 

SA=(SUM-TF)«0,5170 
IF (SA^AHF) 18,17,17 

17 R(N,«JM,0 
GO TO 40 

18 T{N,M)cTF 
TE < N , M ) “SUM 
R(N,M)s3,0 

IF (N-10) 40,37,40 
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37 HMAXrM 
GO TO 40 

30 caul, tlkt,n,m, au,al2,al3,bli#bl2»bl3> 

T(N,M}sT(N,M>»( 1,0*2, 0«AUl^2,0»BLi)*AL2 tt T(N*l l hi>*AU3»T(Nel f M)*BU2« 

lT(N,M+l)*BU3«T(N,M«i> 

40 continue 
50 continue 

IF ITAUwTI) 10.65,65 
C PRINT RESULTS 

65 PRINT 66 i TAU 

66 format <sh time = ,fi0,5,ix,4h min> 

79 FORMAT (UF10.5) 

PRINT 78 

78 FORMAT ( 13H TEMPERATURE ) 

PRINT 72 

72 FORMAT U1H0 N M ,5X,5H TEMP) 

Ns 4 

Ms 5 

PRINT 73,N,M,TlN,M> \ 

73 FORMAT ( 2 I 5 , Fl0 , 5 ) 

N s 5 

M s 9 

PRINT 73,N,M,T(N,M) 

Ns 5 
Ms 13 

PRINT 73,N,M,T(N,M) 

Ns 6 
M s 3 

PRINT 73,N # M,T(N,M) 

Ns 10 

Ms 5 

PRINT 73,N,M,T(N,M) 

NsU 
M s 9 

PRINT 73,N,M,T(N,M) 

NsU 

Ms13 

PRINT 73,N,M,T(N,M) 

Nsl2 
M s 3 

PRINT 73,N,M,T(N,M) 

Nsi6 

Ms5 

PRINT 73,N,M, T(N,M) 

N = 16 

M s 9 

PRINT 73,N#M,f(N»M> 

Ns 17 
Msl3 

PRINT 73,N,M|T(N,H) 

Ns17 
M s 3 

PRINT 73,N,M,T(N,M) 

IF < T I «39 , 9) 82,82,100 
82 T I ?T I +1 ,00 
GO TO 10 
100 STOP 
END 

SUBROUTINE TSHN,M,ASl, AS2 , AS3 , BSl , 8S2 , BS3 ) 

DIMENSION CPS(22,18).AKS(22,18),RH0S(22,18),T<22,18) 

COMMON TO , OX , DY 



Ur'S l PJ | H } s B , 5 l / 

CPS (N^1,M> =0,517 

CP3(N*1»M>=0,517 

CPS (N,M*1>=0,517 

CPS (N,m*1> 80,517 

RHOS(N,M>=0,0160«T(N»M>+54,65 

RHOS(N+1,M>=0,0160*T(N*1»M)+54 I 63 

RHOS(N»*1,M)=0,0160«T(Nb1,K)+54»65 

RHOS(N,M+l>=e,0160«T(N,M*i)+54,65 

RHOS(N»M-i>=0 # 0i60=T(N,M=l)*54,65 

AKS(N,M)=( W ( ,893E^4)»T (N»M) +0,0945) /60,0 

AKS(N*l#M) = (s»( ,893E«4)*T<N+liM>*0,0945>/60,0 

AKS(N-l»M)»(e( ,893E -4 )*T(N-«1#M) *0,0945 )/60,0 

AKS(N,M*l)=<-( ,893E*4 )«TCN»Mr 1>*0,0945)/60 i 0 

AKS(N,M*l)=(»(,893E-4)*T(N,M+l)+0,0945)/60 f 0 

AS1=(TD°AKS(N,M) )/(DX*DX*RHOS(i\i,M)*CPS(N|M) ) 

AS 2 = ( TD*AKS(N*1,M) )/<DX*DX*RHOS(N*i , M ) «CPS ( N* l, H) ) 
AS3=(TD*AKS(N^i,M) )/<DX»DX=RHOS(N-l f M)*CPS(N«l,H) ) 
BSls ( TD*AKS<N,M> >/<DY*DY*RH0S( N,M)*CPS(N,M> ) 
BS2s(TD«AKSCM|H+l) )/(DY*DY*RHOS(N,M*l>*CPS(N|M*l) ) 
BS3s ( TD*AKS<N,M*1) >/(DY*DY*RHOS{N,M«l)*CPS(N»M-i) ) 
RETURN 
END 

SUBROUTINE TS2{N,M,ASl # AS2, AS3#BS1,BS2|BS3) 

DIMENSION CPS(22,18> » AKS(22,18) ,RH0S(22»18> ,T<22|18) 

COMMON TO i OX , DY 

CPS ( N i M > =0 , 517 

CPS(N-l,M)=0,517 

CPS(N+1,K)=0,517 

CPS <N|M*1>*0,517 

CPS(N,M-1) = ,6 057E-3«T{Ni M*1>*0, 46 75 

RHOS ( N, M ) =0 , 0160*T ( N » M ) *54 , 65 

RHOS(N*1,M>=0,0160*T(N*1#M)*54.65 

RHOS CN-l,M)=0,0160ttT(N-l,M >*54,65 

RHOS(N# M* 1>S0,01 60»T(N,M*1>* 54, 65 

RHOS(N,M^i)=«0,0243*T(N,M-J)*50,5 

AKS(N,K )=(»{, 893Ef4)«T(N»M)+0,0945)/60,0 

AKS(N*i,M)s(p ( ,893E«4)*T<N+l#M>*0,0945>/60,0 

AKS(N-l,M>s(»(,893E"A)*T(N^liM)*0 t 0945)/60,0 

AKS(N,M-l)=(-{,893E»A>*T(N#M»l)+0,0945)/60,0 

AKS(N,M*1) = { = ( ,893Ei»A>«T<N#M*l)*0,0945)/60,0 

AS1=(TD»AKS(N,M) )/(DX*DX*RHOS(N,.M)*CPS{N,M) > 

AS28(TD*AKS(N'*1#M) )/(OX#DX*RHOS(N + l,M)«CPS(N*l#M) ) 

AS3=(TD«AKS(N-1,M5 ) / ( DX»I)X*RHOS ( N»l , K ) *CPS ( N»1 » M) > 

BSls( TD*AKS(N,M) )/(DY*DY*RHOS(N,M)»CPS(N,M) > 

BS2={T0«AKS(N,M*1) )/(DY*DY»RHOS(N,M+l)*CPS(N|M*l) ) 

BS38< TD*AKS(N,M«1> )/(DY«DY*RHOS(N,Mttl)=CPS{N»M»l) ) 

RETURN 

END 

SUBROUTINE TU{T,N,M»AU1«AL2»AL3,BL1iBL2»BL3> 
DIMENSION CPU (22, 18) » AKL< 22. 18) , RHOL ( 22 , 18 ) , T ( 22 , 18) 
COMMON TD,DX,DY 

CPU (N,-M)s,6057E»3*T(N,M>+0, 4675 
CPU (N + i,M) = ,6057E«3*T(N*1, M)*0, A 6 75 
CPL ( N-i,M) = , 6057E=.3«T(Nsl, M) +0,4675 
CPU (N#M*1 ) = ,6057E«s3*T(N,K*l)*0, 4675 
CPL(N,Msl) = ,6 357E-*3»T(N#M-i)*0 l 4675 
RHOL(N,M)s W| 0240»T(N»M>*50,5 
RHOl t (N*l,M)s«0,0240»T(N*l l M)*50,5 
RHOL(N»l|M)s*0 l 0248«T(Nsl l M>*50,5 
RHOUIN, M+1)=*0, 0240oT(N, M*l)*50,5 



AKL<N,M)s<»( ,-893E B -4)*T<NiM)*0 t 0945)/60 t 0 
AKL<N*1 ,M)=(b< i 893E«4)»TCN*1,M)*0 i 0945)/60 i 0 
AKL(N-l*M)s(-( ,693E^4)«T(M-1» M> *0,3945 )/60,0 
AKU(N,M*1,) = (»( t 893E-4)«T(N»M*l)+0,09 45)/60 l 0 
AKU(N,Mh 1)={^( t a93Es4)*T(NiM»l)*0 t 0945}/60,0 
ALi?^TD*AKL(N,M) )/<DX»DX*RHQI.(N,M}*CPL(N»M) ) 
AU2 s(T0*AKL(N*i,M) )/(DX*DX*RH0L,(N + 1 i M)»CPL<N*1 i M> > 
Ai-3s(TD*AKL<Nsi l M) >/(OX*PX*RHOL(N»l,M)*CPLCN«.l,Mh 
BU1?CT0»AKL(N.M) )/(DY»DY»RHOU(N,H)«CPU<N»M) ) 
BL2*(TD»AKUN»M + 1> )/(DY*DY*RHOU<N,M*l)«CPL<N»H*l) ) 
BL3 s(TD«AKL(N,M»U )/(DY«0Y«RH0L(N,M*s31)*CPL(N|Ms1> ) 
RETURN 

end 
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CONVECTION MODEL COMPUTER PROGRAM 

This program was written in FORTRAN-IV to solve a 
liquefaction problem of n-octadecane under the influence of 
gravity- induced free convection. The program was run on a 
DEC, model PDP-10, computer. The program was written in 
specific, not general, terms for a given length to width 
ratio and a given number of nodes in spatial directions . 

See nomenclature section for definition of terms. 

Input file: 

1. For execution 4 = Input File or Input Device 

2. First card -TD , DX, DY with a (3F) format where 

TD [=] sec 
DX [=] FT 
DY [=] FT 

3. Second card - ACT, BCT , CCT , DCT , ECT , FCT , with a 
a (6E) format. See 3rd input card section in 
Appendix D for discussion of this input card. 

Flag file: 

1. For execution 6 = Flag File or Flag Device 

2 . Output - 

Time = (time in seconds) 

Input Flag (Hollerith Statement) 

3. Input flag value - (F) format 

a. if flag < 10.0 -continue execution 

b. if flag 2 10.0 -stop execution 

Output File: 

1. For execution 5 = output file or device 

2. Values of velocities for nodal system 

3. time, t, at 120 second intervals 

4. Temperatures of all nodes at 120 second intervals 


Sample Input File: 

First card - 1.0, .005208, .002604 
Second card - -.109E-05, .11859E-03, -.4566E-02, 
. 6941E-01 , - .1157E+00 , .74487E+02 
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NOMENCLATURE 

VMA.X = MAXIMUM VELOCITY ALLOWABLE, FT/SEC 
TD = TIME INCREMENT, SEC 
OX = DELTA X, Ft 

OY = DELTA Y, FT 

• UX S x COMPONENT OF VELOCITY, FT/SEC 

UY = Y COMPONENT OF VELOCITY, FT/SEC 

T = TEMPERATURE AT OL0 TIME STEP, DEG F 
TN = TEMPERATURE AT NEW TIME STEP, DEG F 
TP a TEMPERATURE OF HOT WALL, DEG F 
' to a temperature of cold wall, deg f 
TF = PHASE change TEMPERATURE, DEC F 
AS a SOLID PHASE THERMAL Q I FFUS I V I TY , FT*<»2/SEC 
AL = LIQUID PHASE THERMAL DIFFUSIVITY* FT»»2/SEC 
AHF = LATENT HEAT OF FUSION, R T U / L B 
T I a ELAPSED EXPERIMENTAL TIME, SEC 

r = phase indicator 

+1 IF LIQUID 

-1 IF solid 

DIMENSION U 0 ( 8 2 * 3 4 ) , VC (82. 34 ) 

DIMENSION T (62.34) * T M ( 8 2 » 3 4 ) , R ( 3 2 , 3 4 ) » K ( B 2 1 3 4 ) ,TE(82,34) 
DIMENSION UR( 34 > , NO <82 > 

DIMENSION V V ( 3 4 ) , V ( 82 , 34 ) , U < 82 , 34 ) , V3 < 34 ) , $F < 34 ) , VO ( 34 ) 

C MAXIMUM VELOCITY CALCULATION 
DEL T= 10,0 . 

TD=3,2E-01 

VMAX = 3 ,>224*32. 2* DELT» TO/50. 5 

IF < A8S< VMAX) .GE.2.253E-23) VMAX=2 , 258E-03 

Y 9 5 V M A X 

C COMBINED IDEaL-’VISCOUS VELOCITY CALCULATION 

Y V a - 1 , 0 
8 = 1,0 

As-l,0/SQRT (3.2) 

DO 10 1=2,32 
YV=YV+1.O/16,02 
C = Y V / 8 

RAa(C*tt3.r.C)/<A*»3,-A) 

10 VV < I ) sVMaX»RA 

P I =3 . 1415962 

Y a 0 , 8 

DO 15 1=2,32 

Y = Y + 1 , 0/32 , 88 

X a 1 , 8 - . 

SF ( I ) =S I NH C PI *X > #S I N (PI*Y! 

UBM)a;SINH<PI*X)*C'CS(Pl*Y)*PI 
15 VB< I JacOSH(PI*X)*SlN(PI*Y)*PI 

• X = 0 , 8 

Y = 0 , 58 

DO 18 J=2,16 
X«X*1, 8/32,0 

18 V 0 < 3 4 - j ) = P I * C 0 S H { P I * X ) * S I N < p I * Y ) 

Ya0,50-1. 8/32,8 
DO 36 8=17,32 
Y=Y+1, 0/32,2 
X = 0 , 0 

DO 36 f.=2,16 

X=X*1 , 0/16 , 08 

SFTsSIKH(PI*X)«SIN<PI*Y) 

Ulas-S IKH ( P I *X > *C05 CPI *Y) *.P I 
DO 30 1=1,16 
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KK = 3 3- ! 

IF <SFT,GE,SF(KK) > GO Tq 30 
GO TO 31 

30 CONTINUE 
-KK = 17 

31 UM.AX = UB<KK> 
uaoswtKK) 

U ( N , M ) s Lil » U 8 0 / U M A X 

VN = PI«OCSH(PI»XH>SlN(PI«Y) 

V ( N , M ) s ( 1 , 3 - A 8 S C V N / V B ( K K ) > >*UBO 
IF <ABS<VN'VG(KK> ) .GE.1.00) V<N,M> =0. 0B1*UBO 
36 CONTINUE 

DO 43 h = 2 * 1 6 
DO 40 N s 2 * 1 6 
U(N,M)a-U(N,34-«M) 

V(N,M)*V(N.34*M> 

40 CONTINUE 

DO 5'Z us?, 3? 

DO 50 N s 1 7 i 6 5 

V ( N , M ) = 0 , 1 E - 0 3 
U(N,M)rVV(M) 

IF (M.EO',17) vJ(N.M) s0 t lE-i0 
50 CONTINUE 

DO 60 K‘ a 2 , 3 2 
00 60 M a 6 6 » 8 ‘2 
N N a N - 6 5 

U<81-NN,M)*U(N-64»M) 

V ( 8 1 - N M » M ) 8 • V ( N - 6 4 < M ) 

60 CONTINUE 

C PRINT COMBINED I-V VELOCITY VALUES 
00 62 Ns 2. 80 
DO 62 Mr?, 32 
A U = A 8 S ( U ( N i M ) ) 

A V = A B S { V ( N » M ) ) 

IF ( A U V G E • V 9 ) GO TO 171 
IF UV.GE.V'9) GO TO 171 

62 CONTINUE 

DO 65 M = 2 , 3 ? 

WRITE ( 5 . 63 ) M , ( U { N , M ) * N a 2 » 3 0 ) 

63 FORMAT ( I5i/,5X, ( 6 E) ) 

65 CONTINUE 

DO 70 K a 2 # 3 2 

WRITE (5,681 M » ( v ( N , M ) , N s 2 » 3 0 ) 

68 FORMAT U 5i/#5X,(6E>) 

70 CONTINUE 
GO TO 71 

171 WRITE (6.172) 

172 FORMAT* ' VELOCITY EXCEEDS STABLE VALUE '/) 

GO TO i75 

c read experimental data and input parameters 

71 READ (4,72) TD»DX»DY 

72 FORMAT ( 3F ) 

REAQ (a, 74) ACT, 8 CT,CCT,DCT,ECT,FCT 
74 FORMAT (fiE) 

C SET INITIAL VALUES FOR TEMPERATURE AND FLAGS 

WRITE (5,72) TDiDX.DY 
WRITE (5.74) ACT.BCT,CCT,DCT,ECT,FCT 
T I =0 . 0 £ 

AT 1 = 119,93 
.AHFsl04,90 
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T F s 8 i , 3 3 
T P “ 1 2 0 2 0 
•TO = FCT 
DO 75 m= 2,32 
DO 75 N = 2 i § 2 
NO ( N ) si 
T(N,1)sTP 
T<N,33>=T0 
TN(N.M)=FCT 

TE { Mi N ) =0 . G00 
UC(N#M)s0,C000 
VC(N,M)s0.03 
T(N,M)bTO 

R ( N , M } s - 1 , 0 
K ( N , M > s 1 

75 CONTINUE 

180 T 1 = T ! ♦ T 0 

TA*T 1/60.0 

TO»<<<<ACT*TA+3CT>*TA*CCT)«TA+OCT>*TA*ECT)*TA+FCT 
DO 76 N = 1 1 3 1 

76 T ( M i 3 3 ) = T 0 

C CALCULATE END WALL BOUNDARY VALUES 
J?33 

DO 78 M=2i32 

IF { R ( 2 , M ) , G E , 0 , 1 2 ) GO TO 78 
J = H 

GO TO 79 

78 CONTINUE 

79 DO 30 H s 2 i U 
WsJ,l 

W T = M - 1 

80 T ( 1 * M ) s T p r ( N T / W ) * ( T P - T F ) 

IF ( J , G E » 3 3 ) GO TO 55 

JJ=J+1 

DO 31 m=JJ'32 

W=33-J 

WT*M-J 

81 T<1iM>sTF*<WT/W)*(TF*TO) 

Js33 

85 DO 82 M = 2 # 3 2 

IF <R(fi0,M) .GE.0,10) GO TO 82 

J = M 

GO TO 53 

82 CONTINUE 

83 00 94 Ks2tJ 

Wsj-l 

wt?m-i 

84 T(8l#fi) = TP«(WT/W)*<TP-TF) 

IF (J.GE.33) GO TO 90 

UJs j+i 

DO 86 N= JJ i 32 

W s 3 3 - j 
WTsM-J 
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T < 8 1 , M ) = 

TF* 

-(WT/W)»<TF- 

TO) 

C 

SOLID P 

HA3E 

calculations 
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DO 

112 h 

= 2, 

» 69 



00 

113 M 

= 2 i 

* 32 



IF 

( R ( N i 

M) « 

.GE, 0.1-2) GO 

TO 


SUM 

*0.C 





AS = 

0.975E-06 
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AsT (N,M)e(i t 0-2.3*A3*TD/Qx./QX-2.3*AS«TD/Dy/DY) 
B s T < N * 1 , M ) + T < N - 1 , M ) 

C=AS»T0«3/0X/DX 

DsT (N, m+1)+T(N,M-1) 

E«AS*TO*D/OY/aY 
T N ( N » H ) = A ♦ C ♦ E 

IF <TN(N,M).GE,TF> GO TO 100 
GO TO 110 

C PHASE CHANGE CALCULATIONS 

100 SUMsTN<N,M)-TF 
SUMsSiJM + TE<N,H) 

S A s S U M » 3 , 5 1 7 0 

IF (Sa^aHF) 101.131,132 

101 TE(N,H)=S'JH 
TN(N,H)sTF 
GO TO lig 

102 R(N»M)sl.0 

K < N | M ) s2 

TN(N,K)sTF*(SA*AHF)/0.5l70 

110 CONTINUE 

C LIQUID PHASE CALCULATIONS 
AL*0,950E-06 
DO 115 Ms2,80 
DO 111 1=2,32 
J=I 

IF <R<M, I ) .LT.0,10) GO TO 112 

111 CONTINUE 
fJO(N) »33 
GO TO H5 

112 N 0 ( N ) a j 

115 CONTINUE 

C VELOCIIv CONVERSION CALCULATION'S 
120 DO 132 N; = 2 » 8 3 

L a N 0 ( N ) - 1 
DO 132 1=2,1, 

11 = 1-1 
JJsNQ ( N > -1 

IF (JJ.LT.D JJsMO(N) 

Sl=I I/JJ 

IF ( Si'. GE . 0.999) GO TO 128 
DO 125 J=2 , 32 
JK*J-1 
S2= JK/32 

IF (S2.LT.Si) GO TO 125 
UC(N, I )aU<N, J) 

VC<N, I )=V(N,U> 

GO TO 1,30 
125 CONTINUE 

GO TO 130 

128 UC(N,M)=0. IE-12 

VC<N»M)=(3.1E-12 
130 CONTINUE 

C L.P. CONVECTION ENERGY CALCULATIONS 
DO 135 N = 2 f 3 O 
J s N 0 ( N ) - 1 
IF (J.GE.33) J = 32 
IF (J.LE.2) J=2 
DO 135 H = 2 , J 

IF (K(k,H),GT.1) GO JO 134 
IF ( NO ( N ) . I T , 5 ) U C ( N , M ) s 0 , 1 E - 1 2 



134 

135 
161 

165 
C 

510 

166 

168 

170 

900 

173 

174 


IF ( N 0 ( N ) , L T , 5 ) VC(SiM)=3.1E“12 
A = T (N,M>»( 1,0-2, 3«AL»TD/OX/OX-2.8«AL*TO/OY/OY> 
8 = T (N*l , M ) * < AL»ID/0X/QX-UC < N , M ) *TD/2 . 0/f)X > 

C = T ( Ns-1 , M 5 • ( AL*TD/DX/OX + !JC < N, M ) *TO/P . 3 /OX ) 
D=T(N f M+l)»(AL*TD/DY/DY-VC<N»M)*TD/2.0/DY) 

£ = T AL*TD/0Y/0Y*VC<N.M)*T0/2.3/DY> 
TN(N,M)=A+B+C+G+E 
GO TO 135 
K ( N , M ) b 1 

continue: 

00 165 m=2 , 32 
DO 165 N=2»53 
T(N,M)bTN(N,H) 

PRINT TEMPERATURE PROFILES 
IF < A T I . G C , T I ) GO TO 180 
WRITE (6.510) ATI 

format ( • time = » . F 10 . 2 /) 

A T I s A T I + 1 2 0 , 0 
WRITE (5.166) TI 
FORMAT ( / » F # / ) 

DO 172 Ms 2 , 32 

WRITE (5,168) M , ( T ( N , M ) , N =2,90) 

FORMAT U5,/,<6E)) 

CONTINUE 

WRITE (6,900) <T(N,2),N»l,8l) 

FORMAT ( 4 E ) 

WRITE. (6,173) 

FORMAT (* INPUT FLAG '/) 

READ (6,174) FLAG 
FORMAT (F) 

if (flag. ge, i3,0) go to 175 
IF (TlVLE. 2395,0) GO TO 180 
STOP 
END 


175 


